Criticality in probabilistic models of spreading dynamics in brain networks: Epileptic seizures

The spread of seizures across brain networks is the main impairing factor, often leading to loss-of-consciousness, in people with epilepsy. Despite advances in recording and modeling brain activity, uncovering the nature of seizure spreading dynamics remains an important challenge to understanding and treating pharmacologically resistant epilepsy. To address this challenge, we introduce a new probabilistic model that captures the spreading dynamics in patient-specific complex networks. Network connectivity and interaction time delays between brain areas were estimated from white-matter tractography. The model’s computational tractability allows it to play an important complementary role to more detailed models of seizure dynamics. We illustrate model fitting and predictive performance in the context of patient-specific Epileptor networks. We derive the phase diagram of spread size (order parameter) as a function of brain excitability and global connectivity strength, for different patient-specific networks. Phase diagrams allow the prediction of whether a seizure will spread depending on excitability and connectivity strength. In addition, model simulations predict the temporal order of seizure spread across network nodes. Furthermore, we show that the order parameter can exhibit both discontinuous and continuous (critical) phase transitions as neural excitability and connectivity strength are varied. Existence of a critical point, where response functions and fluctuations in spread size show power-law divergence with respect to control parameters, is supported by mean-field approximations and finite-size scaling analyses. Notably, the critical point separates two distinct regimes of spreading dynamics characterized by unimodal and bimodal spread-size distributions. Our study sheds new light on the nature of phase transitions and fluctuations in seizure spreading dynamics. We expect it to play an important role in the development of closed-loop stimulation approaches for preventing seizure spread in pharmacologically resistant epilepsy. Our findings may also be of interest to related models of spreading dynamics in epidemiology, biology, finance, and statistical physics.


Introduction
Epilepsy is one of the most common neurological disorders affecting approximately 65 million people worldwide [1,2]. About 30% of the cases are diagnosed as pharmacologically resistant epilepsy. Main alternative therapeutic approaches consist of surgical resection of identified epileptogenic brain areas [3,4] or electrical stimulation [5][6][7][8]. In the particular case of focal epilepsy, seizures that initiate in a localized brain region may or may not spread across the brain. While the focal localized seizure onset is commonly not the most impairing aspect, the spreading itself is the main event typically leading to major disruptions in sensorimotor and cognitive processing, as well as loss-of-consciousness. To address the problem of spreading dynamics in epileptic seizures, we introduce a new discrete-state probabilistic network model inspired by continuous-state neural mass Epileptor network models [9]. While complementary to the Epileptor model approach, the probabilistic model allows us here to address fundamental questions about the spreading dynamics which would be otherwise too analytically challenging or computationally intensive.
Data-driven patient-specific Epileptor network models have been widely used to study the dynamics and propagation of focal epileptic seizures [9][10][11][12]. In addition to capturing many dynamical properties of seizures in a small brain area, data-driven Epileptor network models [13][14][15][16][17] have been used to predict seizure propagation in patient-specific networks where each node of the network represents a specific brain area. Dynamics of each node is a neural mass model governed by 6 coupled differential equations (Materials and methods). There are three time scales in the model: a fast time scale for high frequency oscillations during a seizure, a slow time scale capturing both interictal and ictal spike-wave discharges, and a very slow time scale for the variable that captures the dynamics of ionic concentrations and metabolic effects that are thought to be responsible for initiation and termination of seizures. In addition to the above three intrinsic time scales, one can consider an additional refractory period after seizure termination (postictal period).
The connectivity matrices for patient-specific networks are obtained by white-matter brain tractography (structural connectivity) in patients with pharmacologicaly resistant epilepsy. Seizures in the network initiate in specific nodes known as epileptogenic zones (EZs) with high intrinsic excitability where seizures can spontaneously initiate. Beyond its dependence on excitability, seizure initiation in an EZ node depends also on its interactions with non-EZ nodes through diffusive couplings. As a result seizure initiation in the EZ node can be, in some cases, inhibited by non-EZ nodes in a non-seizure state. After initiation in the EZ node, a seizure can spread to the other nodes via the network connectivity, thus affecting the very slow variable of surround target areas through diffusive coupling (Materials and methods).
Depending on the excitability level and connection strengths, the target node may or may not be recruited into the seizure dynamics. These dynamics have been shown to capture the qualitative features of seizure spread in patients with epilepsy [13,14].
Analysis of phase diagrams for seizure spread in patient-specific Epileptor networks, as well as their prediction based on local linear stability analyses, have been examined in [18]. In the control parameter space of excitability and global network coupling strength, three phases are observed: a phase in which the EZ node is inhibited such that not even a seizure in the EZ node is possible (no-seizure), a phase in which a seizure spontaneously initiates in the EZ node but it does not spread (no-spread), and a phase in which the seizure is initiated in the EZ node and it does spread partially or fully (spread) to the surrounding network.
Importantly, as we show in this study, large fluctuations in seizure spread are observed near the transition regions from no-spread to spread phases suggesting critical behavior in the spreading dynamics. However, due to the small size of available patient-specific networks, difficulties in obtaining analytical results, as well as the prohibitive computational cost of simulating very-large Epileptor models on general networks, the nature of these phase transitions and potential critical properties remains unclear. For instance, it is not clear how the size of fluctuations and network responses behave as a function of network size, perturbations in neural excitability and global connectivity strength, or external inputs. The hallmark of criticality, i.e. continuous phase transitions, is the divergence of fluctuation sizes and response functions at the critical point, as well as power-law scaling near the critical point in both equilibrium and non-equilibrium systems [19][20][21][22][23].
As stated earlier, the problem is crucial not only to furthering the understanding of the basic neuroscience of epileptic seizures, but also to the ongoing development of approaches for prevention and control of pharmacologically resistant seizures, such as closed-loop intracranial electrical stimulation via the NeuroPace RNS system. This application context suggests also several choices and constraints to our model development. Given the knowledge that a seizure has just started and remains localized, one would like to determine the optimal intervention (spatiotemporal stimulation) in terms of perturbations of network excitability, connectivity strength, or other related network properties to prevent seizure spreading. Here, we also focus only on typical seizures that self-terminate, i.e. we do not consider the case of status-epilepticus.
To determine the nature of these phase transitions in the above application context, we introduce here a probabilistic trinary network model for seizure spread. Furthermore, we show that the model can be easily fitted to time series data and serve as a computationally efficient alternative for predicting seizure spreading dynamics. Regarding the model fitting, because we find analytic forms for phase diagram boundaries, the model fitting largely reduces to simple curve fitting. Three dynamical states are considered for each network node: susceptible to seizure (preictal), active (seizure, ictal) and postical refractory. State dependent conditional probabilities are defined for transitions between these three states. In contrast with the Epileptor network model, here we are not concerned with the details of the fast dynamics in a node. Instead, we consider only two slow timescales. A time scale for the duration of seizures in a node and a refractory (postictal) timescale. While the three states of the model are reminiscent of common Susceptible-Infected-Recovered-Susceptible (SIRS) epidemiological models [24,25], its dynamics is motivated by that of seizure spread in neuronal networks. In addition, while commonly used SIRS and related models are Markov, the proposed model includes history effects and interaction time delays.
In this study, we focus only on the spreading dynamics within a seizure (our main application context of spread prevention), followed by a refractory or postictal period without transitioning back into susceptible states as in SIRS models. Nevertheless, the model can also be easily adapted to include the inter-seizure dynamics. Because the time intervals between seizures commonly consist of hours if not days, brain areas that have seized will have already roughly recovered to susceptible states (instead of still remaining in the postictal refractory state) before the next seizure. As a consequence, every seizure starts anew in our simulations of the proposed probabilistic model. The initial conditions are always the same: all nodes start in the normal susceptible state, with the surrounding nodes in a non-pathological excitability level, and the EZ area in a hyperexcitable state. Furthermore, parameters such as global excitability, EZ excitability, and coupling strength are assumed to be fixed during a seizure, but are allowed to vary at slower time scales, i.e. across different stochastic seizure realizations. Similarly to nonequilibrium phase transitions in models of directed percolation [20], seizure spreading dynamics in the proposed probabilistic model, under the above specific setup, can have a very large number of absorbing states across different seizure simulations: any spread configuration where a subset of nodes goes into the postictal refractory period.
We show that despite its simplicity, the proposed probabilistic model not only captures seizure spread dynamics and its three phases but also predicts the temporal order of seizure recruitment observed in patient-specific Epileptor network models. In addition, the simplicity and probabilistic nature of the model allow mean-field analysis as well as large-scale computer simulations. Using mean-field approximations and finite-size scaling analyses, we show that transitions from the no-spread to the spread phase include both discontinuous and continuous (critical) non-equilibrium phase transitions, akin to first-and second-order transitions in equilibrium systems, respectively [19]. Furthermore, we show that both the size of fluctuations and response functions diverge at a critical point in the thermodynamic limit. Notably, the size of fluctuations and response functions near the critical point exhibit power-law behavior. At this critical point, the distribution of seizure spread sizes transitions from unimodal to bimodal. These results are supported by numerical simulations of the model on large-scale random networks, and comparisons with Epileptor and proposed probabilistic network models endowed with patient-specific network connectivity.

The model
We model the dynamics of seizure spread in a network of interacting nodes. In the initial following sections, we focus on patient-specific finite-size networks characterized by patient-specific connectivity matrices W, the corresponding interaction time delay matrices τ, and a set of patient-specific epileptogenic and non-epileptogenic zones (EZ and non-EZ nodes) identified by the clinical team involved in the collection of these data (Materials and Methods, and Fig A  and Table A in S1 Text). Specifically, the interaction from node j to i is specified by the interaction weight w × W ij � 0, where w is the global connectivity strength, W ij is the normalized connectivity weight from node j to i, and W ii = 0. Due to the fact that we get the W matrix from white-matter tractography W ij � 0. In addition, we consider time delays τ ij resulting from axonal conduction delays and synaptic dynamics. Each node i is assigned a dynamical trinary variable x i (t) 2 {−1, +1, 0}, where the three values correspond to the three possible states of being susceptible to seizure x i = −1, active (seizure, ictal) x i = 1, and the postictal refractory state x i = 0.
The network system is composed of two types of nodes: EZ and non-EZ nodes. EZ nodes can undergo spontaneous transitions into seizure states, while non-EZ nodes are stable unless their interaction with an active EZ node brings them also into a seizure state. This difference between EZ and non-EZ nodes is reflected in their different neural excitability E j . We assign an E j to each node j depending whether the node is an EZ or not, specifically non À EZ ¼ fj 2 f1; 2; . . . ; NgjE j � 0g: ð2Þ Here we use homogeneous excitability for the non-EZ nodes, i.e. E j = E for all non-EZ nodes. Similarly, for the EZ nodes we set E j = E ez .
In the following, we introduce the stochastic dynamics of the proposed probabilistic network model in discrete time, specified by the corresponding transition rate functions in continuous time. This choice is motivated by our later use of this discrete-time representation in a mean-field approximation of the model's dynamics. Nevertheless, the simulations of the probabilistic network model itself and the direct derivations of its phase diagrams are based on the continuous-time dynamics given by the rate functions. For a small enough time interval Δ, the stochastic dynamics are governed by the following conditional transition probabilities where f and g are rate functions defined below. In words, Eq 3 is the probability of a node i transitioning from susceptible to the active (seizure) state in the time interval (t, t + Δ]. Similarly, Eq 4 is the transition probability from the seizure state to the postical refractory state, and Eq 5 accounts for the transition from the postictal refractory state to the susceptible state.
The rate function f is given by where r is a parameter (in per second). In Eq 3, the argument of f is the total input z i (t) from the network to node i plus the excitability E i of that node. If the input z i (t) to a susceptible node i with x i = −1 is positive and large enough, then there will be a non-zero probability of transition from susceptible to seizure state. Specifically, where the term wW ij u j (t − τ ij ) accounts for the excitatory input from an active node j to a node i that is accumulated over the seizure duration with a time scale τ s . The term u j (t) is formulated so that the accumulated input begins to dissipate with the same time scale after seizure termination. The terms t so, j and t sf, j are respectively the seizure onset and offset times in node j and D j = t sf,j − t so,j is the duration of seizure in node j. The effect of a susceptible node j on node i is conveyed by the term wW ij E j v j (t − τ ij ). Recall the excitability E j is negative for non-EZ nodes and reflects an inhibitory effect. In this way, this term incorporates the diffusive coupling effects featured in the Epileptor model (Materials and methods). The parameters a � 0 and b � 0 account, respectively, for the strength of the two right terms contributing to z i (t) in Eq 7.
In the transition probability from seizure to refractory state, Eq 4, the rate function g is defined as wheret s;i accounts for the time-scale of seizure duration in node i and q s, i specifies the variability of seizure duration. The form of the function g guarantees that seizure duration D i in a node i falls in the ranget s;i À q s;i < D i <t s;i þ q s;i . That is because, if the time since seizure initiation in node i is smaller thant s;i À q s;i , the function gðt À t so;i ;t s;i ; q s;i Þ ¼ 0 (the rate of seizure termination is zero). For a node i, this rate is non-zero in the rangẽ t s;i À q s;i < t À t so;i <t s;i þ q s;i and tends to infinity as t À t so;i !t s;i þ q s;i . Therefore, the duration of a seizure in node i cannot be larger thant s;i þ q s;i . Similarly, for the refractory time T i in a node i we havet r;i À q r;i < T i <t r;i þ q r;i . The variablest s;i , q s,i ,t r;i , q r,i can generally be functions of the history of the system and interactions. Here, we consider the following forms fort s;i and q s,i .t where we assume that the seizure termination probability increases with inhibitory inputs, reducing the expected duration of seizure. The inhibitory effect (cw P which appears in the denominator oft s;i , mimics the phenomenon of surround inhibition in focal seizures which is expected to reduce the seizure duration. This inhibitory effect is assumed to be proportional to the interaction weights wW ij from susceptible nodes which are not in the seizure state. In addition we assume that the variability in seizure duration in each node, specified by q s,i , is proportional to seizure duration. The parameters c � 0 and d � 0 and the function g specify the termination probability of seizures. In Fig O in S1 Text we show that the above choices are in agreement with the seizure duration statistics in patient-specific Epileptor network models (See also Materials and methods: Fitting the proposed probabilistic model to patient-specific Epileptor network simulation data).
As an illustrative example, we show in Fig 1A and 1B that the proposed probabilistic model qualitatively captures the seizure spread dynamics in a patient-specific Epileptor network. In addition, we show how different dynamical variables x i , u i , v i , transition rates gðt À t so;i ;t i ;q i Þ and f(z i + E i ) evolve before during and after the seizure for an individual node in the network. (See Figs B and C in S1 Text for more examples in which the seizure spreads only partially).
In the following sections, we show in more detail that the proposed probabilistic model captures the qualitative nature of seizure spreading dynamics in patient-specific Epileptor networks (Materials and methods). We fit the parameters of the model to capture the same spreading dynamics and corresponding phase diagrams as obtained for the Epileptor network model. We also assume the common scenario seen in actual epileptic seizures where there is no recovery from the refractory state while a seizure is still active in the network. To achieve that, the refractory time scalet r;i ¼ t r is fixed for all the nodes and is much longer than the seizure time scalet s;i .

Phase diagram of the probabilistic model
To derive the phase diagram of the proposed probabilistic model, we considered three main control parameters in the system: the interaction strength w, the excitability of the EZ node E ez and the excitability of non-EZ nodes E. Thus, the parameters (E ez , E, w) did not depend on the network dynamics and were set to fixed values when simulating multiple stochastic realizations corresponding to a given point in the control parameter space for the phase diagrams.
Here we consider a network in which we have a single EZ node ez. The initial condition for the network is {(8i 2 {1, 2,.., N}), x i (0) = −1} where N is the number of nodes in the network. Under this condition, because there is no seizure in the system, all u i (t) = 0 and all v i (t) = 1. Therefore, using Eq 7 we can simply write z ez (t) = Ewb∑ j W ez,j . Next, using the rate function for the transition between the susceptible to the seizure state, Eq 6, it is clear that the transition rate f(z ez (t) + E ez ) is nonzero if Therefore, for the parameter range in the (w, E) space in which this inequality holds, the probability of seizure initiation in the EZ node is nonzero and a seizure will eventually start at the EZ. On the other hand, for the parameter range in which the above inequality does not hold, despite having E ez > 0, the EZ node is inhibited by the susceptible nodes in the network (surround inhibition) that leads to prevention of spontaneous seizures.
Next, we focus on the probability of seizure spread from EZ to non-EZ susceptible nodes. Based again on Eq 6 and assuming that the EZ node is in the seizure state x ez = 1 and that the rest of the network is in the susceptible state x i= 2EZ = −1, we calculated the transition rate of a susceptible node i going to seizure as f(z i (t) + E), where It is obvious that z i (t) + E > 0 leads to a nonzero probability of excitation for node i. Due to the piecewise linear nature of f(y), node i with highest value of z i (t) has the largest probability of excitation. Sitting at the boundary of seizure spread to non-EZ nodes, this node will, most probably, be the first that goes to seizure after the EZ node. Considering the fact that the term u ez (t) is at its maximum when the seizure is about to end in the EZ node, the spread is expected In the above equation D ez is the duration of seizure in the EZ node whose value affects the phase boundary.
The intersection of conditions in Eqs 13 and 15 specifies the parameter range over which there is a nonzero probability of seizure initiation in EZ and a nonzero probability of spread (yellow areas in Fig 2). Additionally, we can specify the parameter range over which there is a nonzero probability of seizure initiation in the EZ node and zero probability of spread (green areas in Fig 2). In this parameter range, the condition in Eq 13 is satisfied, but the condition in Eq 15 is not. Also, there is a parameter range over which the condition in Eq 13 is not satisfied, and regardless of condition in Eq 15, there is no chance of seizure (blue areas in Fig 2). So, we identified three phases of seizure activity in the system: spread, no-spread and no-seizure. Using Eq 13, the boundary between no-seizure and seizure in the EZ node can be evaluated by finding the points in the (w, E) space bellow which the EZ is inhibited and the probability of seizure initiation in EZ is zero (E ez + z ez (t) � 0), and above which the probability of seizure initiation in the EZ node is non zero (E ez + z ez (t) > 0). That can be found by setting (E ez + z ez (t) = 0). Similarly, we can find the boundary between spread and no spread phases using Eq 15. The simulation of the proposed model. Parameters of the model are a = 0.46, b = 0.0021, c = 1.3, d = 0.05, w = 0.45, E = −0.112, E ez = 0.0026. Seizure onset and offset times in each node are shown by a circle and a diamond, respectively. Panels C-F show the evolution of different dynamical variables and transition rates of node 73 (shown in blue in panel B). C Evolution of the state variable x i for i = 73. Before seizure onset x i = −1 (susceptible state); during the seizure x i = 1 (seizure state) and after seizure termination x i = 0 (refractory state). D Evolution of variables u i and v i defined in Eqs 8 and 9. E Evolution of the fraction of active (seizing) nodes in the system and z i (t) as defined in Eq 7. F Evolution of transition rates: f(z i + E i ) for transitioning from the susceptible to the seizure state, and gðt À t so ;t s;i ;q s;i Þ for transitioning from the seizure to the refractory state. Blue dashed lines indicate seizure onset and offset times. For examples of a small spread size see Figs B and C in S1 Text.
https://doi.org/10.1371/journal.pcbi.1010852.g001 boundaries of these regions are respectively for the boundary of no-seizure phase, and for the boundary between spread and no-spread phases. Note that in the above equation i must be the node with maximum z i (t) at t = t so,ez + D ez .

The probabilistic model captures the coarse spreading dynamics in patientspecific Epileptor network models
Here we show that, despite its simplicity, the proposed probabilistic model can capture the seizure spread dynamics of patient-specific Epileptor network models. We first fitted the parameters of the model including (a, b, c, d, τ s ) to time series (summarized as spread size and its temporal evolution) of simulated Epileptor networks. We simulated the Epileptor network model ( Table A in S1 Text for details on patient-specific connectivity matrices and time delays.) In total, we ran 20 different settings. In addition, for each setting, we also varied a range of global coupling w and non-EZ excitability x 0 with 20 stochastic realizations per point in the (w, x 0 ) grid space (Fig 2, top row). We emphasize that we fitted just a single model or set of parameters to the data from all of the simulated patient-specific Epileptor networks and their parameter variations. To fit the parameters τ s , c, d, we examined the seizure duration in the EZ node when there was no seizure spread (green area in Fig 2). In this region of the phase diagram the inhibitory input to the EZ node is constant over time, thust s;ez is independent of time. In this case, we can calculate the probability distribution of seizure duration in the EZ node as a uniform distribution with meant s;ez and standard deviation q s;ez = ffi ffi ffi 3 p (see Materials and methods for more details). Using maximum likelihood estimation, we fitted the parameters c ¼ 1: To fit the remaining parameters in the model, we obtained two different data sets from the simulations of Epileptor network models: (a) the boundary of the no-seizure phase (i.e. the boundary of the blue area in Fig 2), and (b) the boundary between the no-spread and spread phases. Eqs 16 and 17 refer to (a) and (b) above, respectively, and can be used to fit the parameters to the Epileptor network data.
Before fitting the two curves we have also to specify the relation between excitabilities in the Epileptor network model (x 0 , x 0,ez ) and the model (E, E ez ). Considering the excitability level, a single Epileptor network node exhibits a bifurcation point at x 0,b = −2.061. If x 0 > x 0,b , the node is an EZ and can undergo spontaneous seizures; if x 0 < x 0,b , the node is a non-EZ. Our first choice for the relation between the excitability in the proposed model and Epileptor networks is a proportionality (E * x 0 − x 0,b and E ez * x 0, ez − x 0,b ). Here we assume where h > 0 is a parameter to be fitted. Using least-squares, we first fitted Eq 17 obtaining the parameters a = 0.46, b = 0.0021. Finally, we fitted Eq 16 obtaining the parameter h = 0.01 (see Materials and methods for more details).
Using an adaptation of the time-rescaling theorem for point processes [26,27], also known as the temporal Gillespie algorithm [28], we simulated in continuous-time the model with the fitted different parameter sets. The model was simulated with the same patient-specific networks and the same range of parameters (w, E = x 0 − x 0,b ) as for patient-specific Epileptor network simulations. Based on these simulations, we can show that the proposed probabilistic model captures well the seizure-spread phase diagrams of Epileptor networks for different patient-specific network connectivity and different EZ spatial locations (Fig 2 and Fig D in S1 Text).
We note that our simulations of the probabilistic models with the temporal Gillespie algorithm are several orders of magnitude faster than the simulation of the corresponding Epileptor network models. This difference is already clear in the amount of calls to random number generators. While simulation of each realization of an Epileptor network requires two calls to random number generators per node per time step in the numerical solution of corresponding stochastic differential equations, simulation of the proposed model in continuous time via the temporal Gillespie algorithm requires only two calls per state transition in the entire network. For concreteness, the simulation of an 84-node Epileptor network during a 120-second simulation, using a simulation step size of 0.001 s (i.e. M = 120, 000 steps), results in M × N = 1.008 × 10 7 calls to random number generators. In contrast, since a realization of the probabilistic model in continuous time requires only two calls per state transition, the simulation time depends only on the number of transitions between the three states of the model and the number of nodes. For a full seizure spread of size N = 84 we have 2 × N transitions, and thus a much smaller total of 4 × N = 336 calls to random number generators.
Furthermore, the proposed model closely captures the temporal dynamics of seizure spread in patient-specific Epileptor networks (Fig 3). Seizure onset times in the proposed model are linearly related to those of Epileptor networks with a proportionality factor that is different for different points in the phase diagram. As a result of this linear relationship, the order or rank of seizure recruitment in patient-specific Epileptor networks is also well predicted. In addition, the temporal evolution of the fraction of active (seizing) nodes in the network agrees well with the dynamics in the patient-specific Epileptor network models (Fig 4). The above findings generalized over different patient-specific connectivity networks and EZ locations (Figs E-N in S1 Text).
Importantly, as for Epileptor networks, large fluctuations in spread size are also observed for the proposed model near the boundary between no-spread and spread phases ( Fig 5). Studying these fluctuations and their origin is the main goal of the remaining sections of the manuscript.
High fluctuations of this type can appear in near-criticality regions, i.e. near a critical point showing a continuous phase transition. The study of this type of phase transition via numerical simulations is typically computationally very expensive. To overcome this problem, in the next sections we used mean-field approximations to better understand the properties of the phase transition from the no-spread to spread phases.

Random networks and mean-field phase diagram
A mean-field approximation is expected to predict well the dynamics of the probabilistic model on large enough random Erdős-Rényi (ER) networks. In ER networks the connection probability between any two nodes is p. The interaction weight between two nodes is wW ij , where W ij is a random number obtained by sampling from a uniform distribution with mean The plots show the results for patient-specific network P1, EZ-61. A A stochastic realization of an Epileptor network simulation in this patient-specific network. The vertical axis specifies the node index in the patientspecific network while the horizontal axis is time centered at the seizure onset time in the EZ node (node 61). Red diamonds specify the expected seizure onset time predicted by the model. B Linear relation between mean seizure onset times in Epileptor networks versus the mean onset time in the model for all points in phase space in which we observe full spread (yellow area in Fig 2). We note that, while a linear relation between the seizure onset times is observed for all points in phase space, the slope of the line varies depending on w and E. We rescaled all the lines to align with the diagonal. C Seizure-onset ordering in Epileptor networks versus the proposed model for all points in phase space in which we observe full spread. Red dots specify the 95 percentile of the data.
https://doi.org/10.1371/journal.pcbi.1010852.g003 PLOS COMPUTATIONAL BIOLOGY μ 0 /N and standard deviation σ 0 /N. Division by N is necessary for the total interaction strength in the system to be consistent across all network sizes. As a result, the average interaction weight per connection will be wμ 0 /N and the standard deviation wσ 0 /N. Further, for convenience, we set the uniform distribution of W ij to be in the range [0.9, 1.1] × (128/N), such that its mean is 1 for N = 128. Interaction delays are also uniformly distributed in the range [0.75, 1]/60 seconds. In what follows, we set p = 0.2.
We first verify that the mean-field assumption is satisfied for this class of networks in the thermodynamic limit. To be clearer, assume that we have a fraction � of the network as EZ  nodes in the system. The average input weight to a non-EZ node i from all the EZ nodes is The law of large numbers guarantees that the approximation becomes exact in the thermodynamic limit N ! 1 for constant p > 0 and � > 0. As N ! 1, the mean input converges to a constant value wp�μ 0 and its standard deviation converges to zero with � 1= ffi ffi ffi ffi N p . To calculate the phase diagrams in the mean-field approximation, we first generalized Eqs 16 and 17, which were derived for the case where there is only one EZ node in the network, to networks with multiple EZ nodes. For the boundary of no-seizure phase, a generalization of Eq 16 requires that we consider multiple EZ nodes. Since all the EZ nodes have the same level of excitability E ez , the only difference between them is the amount of inhibition that they receive from the susceptible nodes. The EZ node with the smallest inhibitory input is expected to go to seizure prior to the others. Therefore for the boundary of no-spread phase we get where EZg reflects the smallest surround inhibition among all the EZ nodes. Similarly, the boundary between no-spread and spread phases is evaluated under the assumption that all the EZ nodes are in the seizure state and provide excitatory input to other nodes. At the same time each node also receives inhibitory input from the surrounding nodes.
Combining both excitatory and inhibitory inputs, the node with largest z i (t) is the one that specifies the boundary between no-spread and spread. As a result for this boundary we get where the maximization is with respect to time. In the case of just one EZ node (Eq 17), the maximum happens at the time of seizure termination in the EZ node, while in the case of multiple EZ nodes that maximum value can happen at a different time because of different seizure durations. We emphasise that, similar to Eq 17, the index i in this equation indicates the most susceptible node, i.e. the non-EZ node with the highest value of z i (t).
Under the assumption that all the EZ nodes go to seizure state simultaneously, it is straightforward to write the mean-field approximation of Eqs 21 and 22 as and respectively. See Materials and methods for detailed derivations. Eqs 23 and 24 are expected to specify the exact boundaries of the phase diagram for random networks in the thermodynamic limit. While Eq 23 captures well the behavior of finite-size systems, Eq 24 fails to do so (Fig 6). The failure arises from the fact that finite-size ER random networks exhibit, due to variability in degree distribution and also variability in interaction weights, some degree of heterogeneity. The manifestation of this heterogeneity in the dynamics is higher in small networks with a small number of EZ nodes. As a result, the mean-field approximation fails for such small systems.
To address this shortcoming, we provide a finite-size correction that improves the predictions of phase diagrams. This finite-size correction was obtained by multiplying the right hand side of Eq 24 by the correction factor where � denotes the fraction of nodes in the network that belong to the EZ type, p is the random network connection probability, and n denotes a number to be specified of standard deviations in the number of input links from the EZ nodes to a non-EZ node (see Materials and methods for details). As expected, the correction term approaches one as N ! 1. In Fig 6, we show that with n = 2 the finite-size correction captures well the true phase diagrams for different network sizes of N 2 {2 10 , 2 12 , 2 13 }. Red curves indicate the boundary between no-spread and spread phases estimated via mean-field finite-size corrected approximations. Black curves denote the boundary separating the no-seizure phase from the other two phases derived from the mean-field approximations. Bottom row: phase diagrams obtained by the simulation of the mean-field dynamics approximation in discrete-time (Eqs 26, 27 and 28; see also Materials and methods). Red and black curves are the transition boundaries derived from the mean-field approximation without the finite-size correction. As the network size grows, the agreement between the two (top and bottom) phase diagrams improves as expected. https://doi.org/10.1371/journal.pcbi.1010852.g006

Criticallity in seizure spreading: Analysis based on mean-field dynamics
In order to examine the nature of phase transitions and fluctuations in the probabilistic model of seizure spreading, we developed a mean-field approximation of its dynamics. In this meanfield approach, we worked with a discrete-time approximation. Briefly, we express the network dynamics in terms of transition probabilities for the number of nodes entering and exiting different states at specific discrete times. (See detailed derivations in the Materials and methods: Mean-field dynamics section.) In the following, given small enough time intervals Δ, we define the following conditional transition probabilities in discrete time for the non-EZ nodes: PðM mþ1 jH m Þ, Pðn i;mþ1 jH m Þ, and Pðr j;mþ1 jH m Þ. Here, the indices m, i, j denote time bins, and the variables themselves are defined below. The above probabilities can also be similarly defined for the EZ nodes in terms of the corresponding variables M ez i , n ez ij and r ez ij . All these transition probabilities are conditioned on the history of the process up to and including time bin m denoted by where M i and M ez i are respectively the number of non-EZ and the number of EZ nodes that have gone into the seizure state at time bin i. The variables n ij and n ez ij are respectively the number of non-EZ and the number of EZ nodes that have gone into seizure at time bin i and have transitioned into the refractory state at time bin j. Similarly, the variables r ij and r ez ij are respectively the number of non-EZ and the number of EZ nodes that have transitioned into the refractory state at time bin i and have recovered to susceptible state at time bin j.
Specifically, conditioned on H m , the probability of having M m+1 non-EZ nodes transitioning to seizure in the time bin m + 1 can be written in terms of a binomial distribution In the above,t is the mean-field approximation of Eq 11 at time bin m and q s;m ¼ dt s;m is related to the variability of seizure termination times. The term � W denotes the average interaction weight. It is equal to � W ¼ pm 0 =N in the case of ER networks as defined above. Finally, the conditional probability of having r j,m+1 non-EZ nodes (out of all the nodes that have transitioned to the refractory state in time bin j and are still in refractory state) transitioning from the refractory to the susceptible state at time bin m + 1 can be written as wheret r;m ¼ t r is the refractory time scale and q r, m = q r specifies the respective variability in the refractory time. As before, because of our particular setup in this study, r j,m+1 remains zero throughout the simulations. As stated above, the above conditional transition probabilities are also similarly defined for the EZ nodes in terms of the variables M ez i , n ez ij and r ez ij . Details of the numerical simulation of the above mean-field dynamics are given in the Materials and Methods section. We set D ¼ � t=m t , where � t is the average interaction delay computed from patient-specific networks and m τ is the number of simulation time steps in the delay interval � t.
We also derived a correction of the above mean-field dynamics for the case of small networks or sparse spread where the number of susceptible nodes is larger than the number of nodes receiving excitation from the seizing nodes (Materials and methods). That was applied to the case of the small patient-specific networks simulated in this study. We verified that these mean-field dynamics capture well the qualitative features of the temporal evolution of seizure spread size and fraction of active (seizing) nodes as in the probabilistic model and Epileptor networks instantiated in these patient-specific networks (Fig 4 and Fig N in S1 Text). Also, in this case of patient-specific networks, the average connectivity weights � W and the average interaction delays � t, used in the mean-field dynamics approximation (Materials and methods), were computed directly from the patient-specific networks. Specifically, and � t as in Eq 48.
Having derived a mean-field approximation of the spreading dynamics in the proposed model, we proceeded by examining the nature of its phase transitions based on numerical simulations of the mean-field dynamics on ER networks. We examined how the order parameter, i.e. the fraction of spread size s/N behaves as control parameters vary near the boundary between the no-spread and spread phases in the phase diagram. We simulated many stochastic realizations of the mean-field dynamics on a network of size N = 2 15 . Interestingly, we found that the order parameter shows both continuous and discontinuous transitions depending on a specific value of the excitability parameter while varying the global connectivity strength near the phase transition boundary (Fig 7). In particular, there is a clear change of behavior from a continuous to a discontinuous phase transition, suggesting the existence of a critical point. These two different continuous and discontinuous regimes are further manifested in the unimodal and bimodal distributions of seizure spread sizes (Figs 7 and 8). The existence of a transition point between continuous and discontinuous regimes allowed us to approximately locate the critical point.
We examined in more detail the probability distributions of spread sizes via finite-size scaling analyses. The distribution of similar quantities (e.g. avalanche sizes [22]) can show power-law scaling in some, but not all, systems at criticality. Therefore, we checked for the existence of this scaling in the functional form of spread-size distributions near the estimated critical point (Figs 7 and 8). We examined the distributions for increasing network size N near the estimated critical point and also at two other points on the boundary between the no-spread and spread phases. Unimodal distributions, found on the upper part (with respect to the critical point) of the phase boundary, were roughly approximated by Gaussian distributions. Near the estimated critical point, the distributions became more and more skewed, but no evidence of power-law behavior was detected (see also Discussion).
Next, we asked the question of how the size of fluctuations behaves near this critical point. The size of fluctuations at a continuous (critical) transition is known [19,20] to diverge in the thermodynamic limit N ! 1. We consider the standard deviation of spread sizes σ to assess the size of fluctuations in our model. We show that σ appears to diverge at the critical point as a power-law function of excitability and global connectivity strength (Fig 9). The power-law scaling was different on the two sides of the transitions [29]. For fixed E = E c , σ exhibits powerlaw behavior as a function of w as follows: from above the critical point, and from below as The corresponding estimated exponents areĝ ¼ 0:63ð1Þ andĝ 0 ¼ 1:63ð3Þ, respectively (Fig 9A-9C).
Similarly, we examined the behavior of σ with respect to E with fixed w = w c . As before, power-law behavior is also observed for σ according to with estimated exponentsâ ¼ 0:87ð5Þ andâ 0 ¼ 1:4ð1Þ, respectively (Fig 9D-9F). We note that the observed power-law domain increases as N increases, so that fluctuations at the critical point scale as σ m * N 0.66 . Furthermore, the response functions with respect to control parameters are also known to diverge at the critical point in the thermodynamic limit. Here, we defined the response functions of the system as w w ¼ @c @w and w E ¼ @c @E ; where ψ = hs/Ni. Between the red dashed lines we observe bimodality in the probability distribution of the order parameter. The arrow above the critical point indicates a continuous crossover from small spread to large spread sizes. The arrow below the critical point indicates a transition with a discontinuity in the order parameter, i.e. it is not differentiable at that point. Passing through the critical point results in a continuous transition that is expected to exhibit a singularity in the derivative of the order parameter in the thermodynamic limit. We investigated this expected property via finitesize scaling analysis in Figs 9 and 10. Despite the apparent very small region where the above transition from discontinuous to continuous behavior happens, we emphasize that different choices of parameters and their scaling can constrain the seizure spread activity to this small region. For example, based on Eq 21, we note that a choice of smaller EZ excitability (E ez ) level can constrain the spread phase to a very small region around the critical point.
https://doi.org/10.1371/journal.pcbi.1010852.g008 for approaching E c from below, and to w Eþ � ðE À E c Þ 1=dÀ 1 from above. We defined the exponents so that the behavior of the order parameter ψ around the critical point can be written as Exponent symbols were chosen in analogy with the standard use in current literature on criticality. There, the exponent β typically specifies the relation between the order parameter and control parameters, and the exponent δ the relation between the order parameter and an external field. In the standard formulation of critical phenomena, due to the fact that the order Power-law divergence of response functions w w ¼ @c @w and w E ¼ @c @E near the critical point. We used finite-size scaling analysis over four different network sizes of 2 13 , 2 14 , 2 15 , 2 16 . A The expected value of the normalized spread size, ψ = hs/Ni, as a function of w (fixed E) near the critical point (w c � 6.76 10 −5 , E c � 1.00 10 −6 ). The inset zooms the view around the critical point. B The response χ w plotted as a function of w (for fixed E = E c ). The inset shows the divergence of the maximum response χ w,m as a function of χ w,m * N 0.16 (1) . C,D The log-scale plots show the power-law behavior of the response function χ w as w approaches the critical point from below with corresponding scaling χ w− * (w c − w) β 0 −1 and exponent estimated asb 0 ¼ À 0:96ð3Þ, and from above with corresponding scaling χ w+ * (w − w c ) β−1 and exponent estimated asb ¼ 0:43ð5Þ. E The expected value of the normalized spread size, as a function of E near the critical point (for fixed w = w c ). F The response χ E plotted as a function of E (for fixed w = w c ). The inset shows the divergence of the maximum response χ w,m as a function of χ w,m * N 0.18 (1) . G,H The log-scale plots show the powerlaw behavior of the response χ E as E approaches the critical point from below with corresponding scaling χ E− * (E c − E) 1/δ 0 −1 and exponent estimated asd 0 ¼ À 1:6ð5Þ, and from above with corresponding scaling χ E+ * (E − E c ) 1/δ−1 and exponent estimated asd ¼ 12ð5Þ. https://doi.org/10.1371/journal.pcbi.1010852.g010

PLOS COMPUTATIONAL BIOLOGY
parameter is typically zero at the critical point (ψ c = 0 and ψ − = 0), the exponents β 0 and δ 0 are not defined. However, our analysis suggests that these exponents exist in the proposed probabilistic model, where power-law behavior is observed around a nonzero value of ψ = ψ c � 0.14.
In other words, the estimated critical point lies inside the region of the defined spread phase. This can happen due to finite-size effects in some models, e.g. kinetic Ising models [30] and epidemic models in the dynamic isotropic percolation universality class [31]. In such models, the estimated ψ c will decrease and approach zero by increasing system size. Nevertheless, it seems that this is not the case in the proposed probabilistic model because we observe that ψ c does not decrease by increasing system size. We also note that the bimodal regime occupies a thin region almost parallel to the phase boundary (Fig 8I and 8J). The unimodal regime, with the mode consisting of either a small or large spread, is observed everywhere else inside the spread phase region.
The numeric values of these exponents were estimated asb 0 ¼ À 0:96ð3Þ,b ¼ 0:43ð5Þ, d 0 ¼ À 1:6ð5Þ, andd ¼ 12ð5Þ. While the exponentsb andd are positive and in the expected range, the estimated value ofb 0 andd 0 are negative. This negative exponent leads to a singularity in ψ and that is not acceptable as 0 � ψ � 1. We think that these inconsistent exponents result from what has been referred before as apparent exponents [32]. They appear when the scaling function exhibits power-law behavior in such a way that masks the actual critical exponent. (See Materials and methods for a formulation of the behavior of ψ − in terms of apparent exponents.) All of the above estimated exponents are summarized in Table 1. In sum, to our knowledge, the proposed probabilistic model does not belong to any of the well known universality classes. We note, nevertheless, that the estimation of critical points is prone to finite-size effects and numerical inaccuracy. It is possible that more accurate methods may lead to slightly different exponents. We hope these estimated exponents will help to shed some light on the spreading dynamics of epileptic seizures. In particular, the exponents inform about the sensitivity of the modeled spreading dynamics to perturbations in the control parameters (w, E) and also in external inputs to the system near the critical point.
As a complementary evidence of criticality, in addition to the above described power-law behavior of response functions, we also found that the maximum values of response functions χ w,m , χ E,m tend to diverge with increasing network size N. Furthermore, these values also followed power-law functions of N according to χ w,m *N 0.16 (1) and χ E,m * N 0.18 (1) . Table 1. Estimated exponents. The terms including σ denote the standard deviation (size of fluctuations) of spread size and their dependence on excitability (E) and connectivity strength (w) in the probabilistic model. The terms including ψ relate to the order parameter (normalized spread size) and their dependence on excitability and connectivity strength. The numbers in parentheses after each value represent the error in the last digit, e.g 1.63(3) = 1.63 ± 0.03.

Discussion
The nature of the spreading dynamics in epileptic seizures remains a challenging problem in neuroscience, with important implications to the development of new therapeutic approaches, especially in the case of pharmacologically resistant seizures. Here, we have provided two novel contributions to address this challenge. First, we have introduced a new probabilistic network model that can capture the spatiotemporal spreading dynamics of seizures in patientspecific complex brain networks. We set up the problem in the context of a focal seizure that has just started and initially remains localized into the epileptogenic zone. The question is then whether the seizure will spread or not, and how. To answer this question, we started by showing that the model can be fitted to data-driven patient-specific Epileptor networks. Because of its probabilistic and phenomenological nature, the model can be easily fitted and is fast to simulate. We then derived the phase diagrams of patient-specific models, where the order parameter was defined as the seizure spread size, and the control parameters were defined as the neural excitability and global connectivity strength. The phase diagrams allowed us to determine whether a seizure will spread based on the excitability and global connectivity strength in the brain. We have also shown that simulations of the model also successfully predicted the temporal evolution of spread size and the temporal ordering or rank of different network nodes in the surrounding as they are recruited into seizure. In this way, fast simulations of the probabilistic model can accurately predict how a seizure spreads. Second, our analyses revealed the nature of the phase transitions in the seizure spreading dynamics in these probabilistic models. We have demonstrated that the order parameter spread size can show both discontinuous and continuous (critical) phase transitions as neural excitability and global connectivity strength are varied. A mean-field approximation of the dynamics and finite-size scaling analyses provided supporting evidence for the existence of a critical point near the boundary separating the no-spread and spread phases. Specifically, we have shown that the standard deviation of fluctuations in spread size diverges with power-law scaling at numerically estimated critical points as the network size N is increased. Furthermore, we have also shown that the corresponding response functions, namely the partial derivatives of the order parameter with respect to excitability and global connectivity strength, also diverge at the critical point with increasing N. Importantly, this critical point separates two distinct regimes in the spreading dynamics for control parameters near the boundary between no-spread and spread phases. These two regimes are characterized by either unimodal or bimodal probability distributions of spread size. In the unimodal regime, seizures present spread sizes that range from very small to large number of network nodes in a continuous fashion. On the other hand, in the bimodal regime, seizure spread sizes are typically either very small or very large, with few or no observed intermediate sizes. Stochastic fluctuations trigger either small or these seemingly explosive large spread size events.
We emphasize that currently there is not enough data from recordings in either patients or animal models, both in terms of the number of recorded seizures per subject and in terms of recordings with full coverage of brain areas in both hemispheres. The number of recorded seizures per subject tends to be very small, especially in the hospital setting of epilepsy monitoring units. Available ECoG or SEEG recordings have commonly been restricted to a small subset of brain areas candidate for resective surgery or device implantation. All of this makes the examination of distributions of seizure spread sizes and related statistics based on experimental datasets currently unfeasible or at least very incomplete. Given these data restrictions, here we have adopted patient-specific Epileptor networks as our main reference for the development of the probabilistic model of seizure spreading. As stated earlier, patient-specific Epileptor networks have been fitted to actual patient data and shown to successfully capture many of the seizure dynamics features [13,14,17,33]. In this sense, we think patient-specific Epileptor networks provide an initial good reference for the study of probabilistic models of spreading dynamics in epileptic seizures. Nevertheless, we also note that, given the complexity of Epileptor network models and their costly computational simulations required in finite-size scaling analyses, for example, more concise probabilistic models as the one proposed here can play a fundamental complementary role. We should also emphasize that since the proposed model is a coarsened version of the Epileptor network states, it does not capture all their dynamical features such as synchronization of oscillatory activity across networks nodes and its potential effects on seizure spread. In this sense, we do not claim our analyses elucidate the nature of the phase transitions in the Epileptor model in its entirety. That remains an open question. The presented mean-field analysis is a mean-field analysis of the proposed probabilistic model.
As stated above, our focus in this study has been the scenario more immediately relevant to predicting whether and how a seizure that has just started will spread. In this first treatment, we have ignored how inter-seizure dynamics and potential history effects of spreading patterns in a given past seizure might affect the spreading in future seizures. In other words, we assumed the spreading dynamics in a given seizure to depend only on the within seizure history related to the events that have happened since the seizure initiation and that are independent from spreading events in previous seizures. This might be a reasonable first-order approximation if seizures are sufficiently far apart in time. In terms of model improvement, inter-seizure dynamics can be incorporated via the inclusion of appropriate time scales for the recovery from postictal periods. In addition, recent studies suggest that the temporal dynamics within sequences of seizures show multiple time scales, from circadian to multiday rhythms [34]. Furthermore, once these dynamics are incorporated, one would also need to address the possibility of seizure spread patterns in a given seizure to affect the spreading pattern in future seizures. This might involve synaptic plasticity and maintenance of pathological networks for seizure spread.
Given the choice of parameters and existence of fixed epileptogenic areas that start seizures and drive spreading, and the specific allowed sequences of events to capture typical neural dynamics during seizure onset and spread, the proposed probabilistic model lacks detailed balance, the condition for thermodynamic equilibrium. Therefore, our use of phase transition and related terminology needs to be considered outside the statistical physics of systems in equilibrium. The concept of phase transitions and critical behavior has over the years been extended to systems near and far from equilibrium, e.g. [20,35]. The main related reference to our study here is the field of directed percolation (DP) and other related models of nonequilibrium systems [20,[36][37][38][39]. Although the proposed probabilistic model is more complex and does not belong to the DP universality class, our use of concepts such as phase transitions, criticality, and our power-law scaling analyses are analogous to that commonly performed in the DP and related fields.
We have demonstrated the existence of signatures of criticality based on divergence of fluctuations in spread size and related response functions at estimated critical points. Furthermore, the exponents for power-law scaling differed between the two sides of the transition [29]. However, in contrast to many previous studies focused on distributions of avalanche sizes, our analysis did not find power-law behavior in the functional form of distributions of spread size at the estimated critical points. Nevertheless, we emphasize that, although this is expected in many processes, e.g. avalanches in self-organized critical systems [21,22], powerlaw scaling in the probability density of the order parameter (here spread size) can also be absent in many other systems at criticality. Examples of the latter can be found in Ising models where the net magnetization can show bimodal distributions with no power-law scaling [40,41].
The estimated critical point inside the spread phase region, where fluctuations and response functions appear to diverge, is also at the transition between the two regimes of unimodal and bimodal distributions of spread size in the examined ER random networks. Two main related issues can be raised. First, this transition suggests a bifurcation from a single (small spread) to two metastable states (small and large spread) in the spreading dynamics. One could argue that our finite-size scaling results would be more properly interpreted from the perspective of bifurcation analyses in dynamical systems under stochastic perturbations, rather than from the statistical physics perspective of phase transitions. Nevertheless, these two perspectives are not necessarily mutually exclusive and can be applied in a complementary manner when approaching some nonequilibrium systems, e.g. [35,42]. Second, it has been shown that in mean-field analyses of certain stochastic nonlinear systems, e.g. stochastic Wilson-Cowan dynamics on random networks, metastability vanishes in the thermodynamic limit N ! 1, being replaced by multistability in the resulting deterministic system [43]. In addition, the transition rates between metastable states are expected in this case to decay exponentially with increasing N. We did not observe such behavior in our simulations. In contrast, bimodality was enhanced with increasing N. Furthermore, we emphasize that in this bimodal regime where stochastic realizations can result in either small or large spread, the occurrence of the later always requires the network to approach first the small spread metastable state. That results, in our setup, from the initial conditions of the network always being the normal susceptible state, i.e. zero spread size. Exponential decay of transition rates would imply that the occurrence of large spread sizes should become less and less likely with increasing N, something that we also did not observe. Another potential issue that could be raised is that this bimodal regime could somewhat contribute artifacts to our finite-size power-law scaling analyses. We think that such effects are unlikely given that our finite-size analyses are performed at the estimated critical point and for directions including only the unimodal regime. In sum, although these issues remain open to further investigation, we think that the above arguments support the approach taken here.
While our initial results relating the probabilistic model to Epileptor networks were based on patient-specific complex networks, i.e. networks with modular graphs reflecting the hemispheric and other brain areas organization structures, our mean-field approximations and results relied on the assumption of Erdős-Rényi random networks. Although this disparity did not affect significantly the prediction of phase diagrams and the time evolution of seizure spread size in patient-specific Epileptor networks (Figs 4 and 6, and Fig N in S1 Text), the extension of our response function and finite-size scaling results based on the mean-field approximations to these complex networks remains an important open problem. It also remains to be shown how our results extend to patient-specific networks with much finer parcellation of brain areas and more varied epileptogenic zone locations than those examined in this study. We hope to address these problems in future studies.
We expect that the type of probabilistic models for spreading dynamics proposed here will play a fundamental role in closed-loop intracranial stimulation control approaches for preventing seizure spread. For instance, predictions based on these patient-specific probabilistic models can guide the specification of spatiotemporal stimulation patterns in NeuroPace RNS devices [7,44] endowed with the ability to track biomarkers of brain excitability and connectivity strength. Finally, we hope that the continuing development of new human neurophysiological recording strategies and devices will soon allow for the experimental testing of the many predictions derived by our analysis of seizure spreading dynamics.

The Epileptor network model
We follow closely the formulation in [15,18], keeping the notation for the Epileptor network model the same as in previous publications so that they can be easily related. Some of the symbols overlap with the notation for the proposed probabilistic model, but the distinction should be clear from context.
For an N-node patient-specific Epileptor network model, the dynamics at each node i = 1, 2, . . ., N are given by where gðx 1;i Þ ¼ The terms ξ i (t) and η i (t) correspond to zero mean Gaussian white noise. Also, both ξ i (t) and η i (t) are independent across nodes. Here, we set the corresponding noise variances to 0.0025. Effectively, we interpret the above as stochastic differential equations in the Itô calculus sense and apply the Euler-Heun method in the stochastic simulations. We set I 1 = 3.1, I 2 = 0.45, γ = 0.01, τ 0 = 6667, τ 1 = 1, τ 2 = 10. For agreement with previously published simulations of Epileptor network models, we list the time related variables in time units of 0.02 seconds, e.g. τ 1 = 0.02s.
The coupling weights W ij � 0 were obtained from patient-specific connectivity matrices derived from white-matter tractography. The terms τ ij > 0 are the corresponding axonal transmission delays. (For more details, see Structural network connectivity below.) The global connectivity strength is specified by the parameter w. We note that, although the connectivity matrix and global coupling consist of nonnegative or strictly positive values, network interactions can result in suppressive or inhibitory effects via a diffusive coupling between nodes instantiated as [x 1,j (t − τ ij )−x 1,i (t)] in the above equation for the z variable.
Epileptogenic and non-epileptogenic nodes can be instantiated by setting the corresponding node excitability parameter x 0,i to specific distinct values. This excitability parameters control whether a node can go spontaneously into seizure. Its critical value for an isolated node is x 0,b � −2.061. Levels above this value allow for spontaneous seizure transitions. Here, for an epileptogenic zone (EZ) we set x 0,i = {−1.6, −1.8}. For non-epileptogenic nodes in the surround network we considered x 0 2 [−2.3, −2.09], with 0.0115 steps.
For a detailed motivation of the Epileptor network equations and their dynamics we refer to [9,10,33,45,46]. Briefly, three different time scales, reflected in the parameters τ 0 � τ 2 � τ 1 , allow the network to capture both the slow and fast oscillations typically observed in epileptic focal seizures. The slowest time scale corresponds to potential ionic, metabolic and homeostatic dynamical processes leading the network to transition into seizures. This slow dynamics is captured by the "permitivity" variable z i (t), which effectively works as a bifurcation parameter, leading the network to spontaneously transition in and out of seizures.

Patient-specific structural network connectivity and interaction delay matrices based on white-matter tractography
Details about (diffusion MRI) white-matter tractography, brain area parcellation, and plots for all of the 5-patient network connectivity and time delay matrices used in this study are provided in Fig A and Table A in S1 Text.

Fitting the proposed probabilistic model to patient-specific Epileptor network simulation data
The proposed model has in total 6 parameters {a, b, c, d, τ s , h}, which need to be estimated. To fit these parameters to data from patient-specific Epileptor network models, we simulated Epileptor networks instantiated on 5 different patient-specific networks. For each network, we considered two different cases of EZ node location. For each network, we also considered two different EZ excitability levels set at x 0,ez = {−1.6, −1.8}. In total we had 20 different networks. In addition, we varied the control parameters (w, x 0 ), i.e. the excitability and global connectivity strength, respectively, on a grid space. For each point on these grid spaces and for each network, we simulated 20 stochastic realizations and computed the corresponding phasediagrams as shown in Fig 2. As stated earlier, we emphasize that we fitted just a single model to all of the simulated patient-specific Epileptor networks, choices of EZ node location, and variations in excitability and global connectivity strength. First, to fit the parameters {τ s , c, d}, which specify the duration of seizures and their variability, we considered the seizure duration in the EZ node D ez in the no-spread phase region of the (w, E) space. In that region, the seizure starts in the EZ node but does not spread (green area in Fig 2). Furthermore, in this region of the phase diagram, the input to the EZ node is constant over time Thus, we can simply calculate the probability distribution function of D ez during this phase of the network using Eqs 5, 10, 11 and 12 by considering y = t − t so,ez as the time elapsed since the seizure initiation, the relation between the duration of seizure and rate of termination is pðD ez Þ ¼ gðD ez ;t s;ez ; q s;ez ÞexpðÀ R D ez 0 gðy;t s;ez ; q s;ez ÞdyÞ: Because in this parameter ranget s;ez and q s,ez are independent of time, direct calculation leads to a uniform distribution of D ez in the range ½t s;ez À q s;ez ;t s;ez þ q s;ez Þ with meant s;ez and standard deviation s D ez ¼ q s;ez ffi ffi ffi 3 p : Next, we examined the behavior of D ez as a function of the input to the EZ node (Fig O  panel A in S1 Text). Furthermore, the expected value conditioned on I ez and standard deviation of D ez are both in agreement with the functional form of Eqs 11 and 12 (Fig O panels B and C in S1 Text). By fitting the functions directly to the data, we estimated τ s = 32.22s, c = 1.3 and d = 0.05.
Having determined τ s , to fit the parameters a, b we extracted the phase transition boundary (the boundary between green and yellow in Fig 2) from all of the phase diagrams constructed for all patient-specific networks and corresponding EZ nodes from the Epileptor simulations. For each phase diagram, the data were constructed as the set of pairs {(E k,r , w k,r )}, where the indexes k and r indicate different points on the given transition boundary and different stochastic realizations, respectively. We used least-squares to fit parameters a and b. Specifically, we used Eq 17 in the following quadratic cost function to be minimized with respect to a and b leading to a = 0.46 and b = 0.0021. We note that the index i in the above cost function specifies the most susceptible node, which is identified by finding the node with maximum input from the EZ node. In the above equation the first and second summations are with respect to patients and EZ nodes, respectively. In order to fit the remaining parameter h, we extracted the "inhibitory" boundary (i.e. boundary of the blue region in Fig 2) from each of phase diagrams obtained from the Epileptor network simulations. For each case, the data were constructed as the set of pairs of {(w j,r , E j,r )} where the indexes j and r indicate different points on the phase boundary and different stochastic realizations, respectively. Using Eq 16 and setting we minimized the following quadratic cost function resulting in h = 0.01.

Numerical simulations: Temporal Gillespie algorithm
To simulate in continuous-time the exact dynamics of the probabilistic model we used an algorithm based on an extension of the time-rescaling theorem [26,27], also known as the temporal Gillespie algorithm [28]. Consider first just a single-node (i.e. N = 1) with a corresponding set O of possible state transitions including seizure initiation, termination (transition to refractory or postictal) and recovery. Any realization of the process consists of a sequence of transition times ft 1 ; t 2 ; t 3 ; . . .g; which can be represented as a realization from a stochastic point process with conditional intensity where HðtÞ is the history of the process, and l o ðtjHðtÞÞ is the conditional intensity for each of the possible ω 2 O transition types. The l o ðtjHðtÞÞ are specified by the corresponding rate functions f and g, i.e. Eqs 6 and 10, with appropriate input arguments. By the time-rescaling theorem, the rescaled (waiting) time intervals between transitions are independent and exponentially distributed with mean 1, i.e. the rescaled-time point process corresponds to a homogeneous Poisson process with unit rate. Thus, to simulate the original process, one can simply start by sampling a sequence of waiting time intervals {υ k } from this unit mean exponential distribution, and for each one of them, use Eq 39 to solve for the transition times. Concretely, given a transition time t 1 and a new sampled waiting time υ 1 , one solves Eq 39 for the next transition time t 2 . Next, the type of the transition itself, indicated by ω, is sampled according to the probabilities The simulation continues from one transition to the other by repeating the above steps. The above can be easily extended to networks, by setting where λ j,ω corresponds to the intensity at node j, and using again Eq 39 to solve for the transition times. Finally, both the node and the type of the transition are sampled with probabilities where t k is the currently sampled transition time under consideration. The simulation proceeds following the above steps. We emphasize that the seizure always starts in the EZ nodes. Furthermore, when applying the algorithm to ER random networks (as opposed the patient-specific networks), we enforce a seizure to start simultaneously in all of the N ez = �N EZ nodes.

Calculation of the spread to no-spread boundary in the mean-field approximations
To evaluate the boundary between seizure spread and no-spread, we start with the case where the EZ nodes are in the seizure state. We want to find the points in the parameter space (E, w) at which the probability of spread is zero. We refer to the boundary of this region as the no-spread to spread boundary. Assuming that the EZ nodes are in the seizure state, following the same steps as in deriving Eq 17, and considering the node with largest z i (t), we can write the equation for this boundary as While the above equation is written for the node with largest z i (t), in the mean-field approximation all the surrounding nodes receive equal input from the EZ nodes. Therefore we are not concerned with the maximization with respect to i. To derive the mean-field approximation of the boundary between no-spread and spread phases, we replace W ij with its average value � W and replace the interaction delays τ ij with their average value � t. Additionally, since the the summation in the denominator of the above equation does not include the EZ nodes, it can be evaluated as P Next, to get the final form of this boundary we need to find the maximum of with respect to time. Note that as the translation in time with � t does not change the maximum value of this function, for convenience we do not consider it in the following calculations. In our setup, we consider the case that all the EZ nodes go to seizure simultaneously at time zero.
(This reflects the observation that the time scale for seizure spread within the EZ region is much faster than the time scale for seizure spread in the large-scale network.) In addition, in the no-spread phase, the duration of seizures in the EZ nodes is a random variable uniformly distributed in the range ðt s;ez À q s;ez ;t s;ez þ q s;ez Þ. Based on the above, we can evaluate an approximation of U(t) that is exact in the N ez ! 1 limit. To achieve that, we divide the time axis into 4 intervals, and for each interval we calculate the function U(t). See Fig P in S1 Text for visualization and a summary of the following results. The intervals and related calculations are: 1. Time interval 0 � t <t s;ez À q s;ez : All the EZ nodes are in the seizure state and all have the same u i (t) value. Therefore, 2. Time intervalt s;ez À q s;ez � t <t s;ez þ q s;ez : The uniformly-distributed (for large N ez ) seizure termination times for EZ nodes can be approximated as homogeneous termination times, which are proportional to the time passed sincet s;ez À q s;ez . In other words if �ðtÞ ¼ t À ðt s;ez À q s;ez Þ and F = 2q s,ez , then the ratio of EZ nodes that have undergone seizure termination up to time t is equal to �ðtÞ F , and the ratio of nodes that are still in the seizure phase is 1 À �ðtÞ F . In this time interval, the contribution of the nodes that are still in the seizure state to U(t) is t t s 1 À �ðtÞ F À � . Using the fact that for the rest of the nodes the {u i (t)} are in their linear decreasing phase, and also the uniform termination times, the average contribution of these nodes is equal to �ðtÞ F ðt s;ez À q s;ez Þ. Therefore, Since for t >t s;ez þ q s;ez all fu i ðtÞji 2 EZg are in their decreasing phase, the maximum point of U(t) is expected to be in the above intervals. In addition, the maximum point is not in the first interval because UðtÞ < Uðt s;ez À q s;ez Þ for any t <t s;ez À q s;ez . Therefore, the maximum point must be in the second interval. By calculating the root of the time derivative of U(t), which is a quadratic function with a unique extremum point in the interval, we can show that the maximum point happens at t ¼t s;ez . Thus U max ¼t s;ez À q s;ez =2 t s : For completeness, to show the calculation of U(t) at all times, we also evaluate it here in the following two intervals.
3. Time intervalt s;ez þ q s;ez � t < 2ðt s;ez À q s;ez Þ: The {u i (t)} for all the nodes are in their decreasing phase, and due to the uniform distribution of the termination times, their average is equal to 4. Time interval 2ðt s;ez À q s;ez Þ � t < 2ðt s;ez þ q s;ez Þ: The {u i (t)} of some of the nodes are in their linearly decreasing phase and u i (t) = 0 for the rest of the nodes. The ratio of the nodes for which u i (t) = 0, during times in this interval, is proportional to the time that has passed since 2ðt s;ez À q s;ez Þ, i.e. tÀ 2ðt s;ez À q s;ez Þ 2D . The ratio of the nodes that are in the decreasing phase is 1 À tÀ 2ðt s;ez À q s;ez Þ 2D . Therefore, the average u i (t) of these nodes is equal to UðtÞ ¼ ð1 À t À 2ðt s;ez À q s;ez Þ 2F Þ 2ðt s;ez þ q s;ez Þ À t 2t s : Outside of the union of these four intervals, u i (t) = 0 for all the EZ nodes and obviously U(t) = 0.
Having derived the maximum value of U(t) we continue the calculation of the boundary of spread to no-spread. By plugging U max ¼ max t f P j2EZ u j ðtÞ=N ez g ¼ ðt s;ez À q s;ez =2Þ=t s into Eq 41, we get the boundary of no-spread to spread as Based on the mean-field approximation for large random networks with connection probability p, with a fraction of EZ nodes � = N ez /N and connectivity weights sampled from a distribution with average μ 0 /N, we get � Finite-size correction for mean-field derived phase diagrams We present the detailed derivation of the finite-size correction Eq 25 for the mean-field prediction of phase diagrams for the proposed model. Having a fraction � of all the nodes in a random network, with connection probability p, as EZ nodes, the total number of edges from all the EZ nodes to a node in the surrounding is distributed according to a binomial distribution with mean Multiplying the total number of edges by the average connection weight per edge (μ 0 /N), we get the average and standard deviation of the input from EZ nodes to a random node in the network as respectively. For finite N there is variability in the inputs from the EZ node to surrounding nodes, and there is a probability of having non-EZ nodes that receive inputs larger than � W EZ . A finite-size correction to Eq 23 can thus be obtained by considering such nodes as the most susceptible nodes with inputs up to n standard deviations larger than the mean input, leading to Therefore the finite-size correction is obtained by multiplying the right hand side of Eq 24 by the correction factor n ¼ 1 þ n ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ð1 À �Þ �pN s :

Mean-field dynamics
Mean-field analysis is based on the assumption that the effective input to a node from others at a given time is the same for all the nodes in the system. In other words, if we define then under the mean-field assumption � W i ¼ � W j for all pairs of the nodes. Therefore, we drop the index of this mean value and use � W . Now, we can write the mean-field expression for z i (t) as where � t is the average interaction delay over all pairs of nodes defined as where A ij are binary elements of the adjacency matrix (i.e. A ij = 1 if there is a connection form j to i, and zero otherwise). Using the above equation and Eqs 3-5, we derive the mean-field equations for evolution of the number of nodes in each of the three possible states of the model.
To do that, we make a discrete-time approximation and study the evolution of the system over small time intervals of size D ¼ � t=m t , where m τ is the number of time steps in the delay time interval. We define the number of non-EZ nodes that transition to the seizure state at each time bin i as M i , the number of non-EZ nodes that enter the seizure state at time bin i and exit this state at time bin j (seizure ending) as n ij , and the number of non-EZ nodes that exit the seizure state in time bin i and recover from refractory state in time bin j as r ij . Similarly, we define the same variables for the EZ nodes as M ez i , n ez ij and r ez ij . Next, using Eqs 8 and 9, we can write z m = z(mΔ) in terms of the above variables as in Eq 49 gives the number of nodes that have entered the seizure state at the time bin i and are still in this state at time bin m. It accounts for the linear increase of the output from these nodes to the neighbor nodes before seizure termination. On the other hand, the term X m j¼i Hð2j À i À mÞðn ij þ n ez ij Þ accounts for the linear decrease in the accumulated input after seizure termination (see Eq 8).

The function
HðxÞ ¼ We write the above probability function based on the fact that there are N s, m susceptible nodes at time bin m, out of which M m+1 nodes are chosen to transition to the seizure state at time bin m + 1 where each transition takes place with probability f(z m + E)Δ.
Similarly, we can write the probability of having n i, m+1 non-EZ nodes, from the M i nodes that have gone to the seizure state in the time bin i and are still in the seizure state, transitioning to the refractory state at time bin m + 1 as is the mean-field approximation of Eq 11 at time bin m and q s;m ¼ dt s;m is related to the variability of seizure termination times. In the above equation � W is the average interaction weight. Finally, the probability of having r j, m+1 non-EZ nodes, out of all the nodes that have transitioned to the refractory state in time bin j and have not recovered yet, transitioning from the refractory to the susceptible state at time bin m + 1 can be written as Pðr j;mþ1 jH m Þ ¼ respectively. Eqs 51-53 can be similarly written for the EZ nodes represented by the variables M ez i , n ez ij and r ez ij defined above.

Sparse-seizure correction to mean-field dynamics
We developed a correction to Eqs 26 and 49 that is valid in the regime of spreading dynamics where there is a small number of output edges from seizing nodes to susceptible nodes compared with the total number of susceptible nodes in the network. Examples of this regime are the early time of spread dynamics on networks with small number of EZ nodes, e.g. patient-specific networks examined here, and sparse random networks with small average degree. For a random network with connectivity probability p, the spread dynamics at time bin m is determined by the number of seizing nodes N e,m and the number of susceptible nodes N s,m . The expected number of output links projecting from the N e,m active nodes to N s,m susceptible nodes can be written as K se,m = pN e,m N s,m , where the subscript se indicates the links or edges from excited to susceptible nodes. If N s,m > K se,m , then the output from the active nodes does not evenly distribute among all the N s nodes; instead it can at most affect K se,m of the N s,m susceptible nodes. Therefore, we can make a correction to Eq  ! .
In addition, we can also make a correction to Eq 49. Based on the above argument, the excitatory output from the N e,m nodes can at most excite K se,M nodes. As a consequence, we can rescale the excitatory part of z m in Eq 49, which is evaluated to distribute among all susceptible nodes by a factor of N s;m K se;M : Obviously, the above corrections are valid only if N s,m > K se,m . Otherwise, the dynamics follows the originally derived mean-field dynamics. We used mean-field simulations with the above corrections to approximate the spreading dynamics in patient-specific networks as shown in Fig 4C.

Simulation of mean-field dynamics
To simulate the mean-field dynamics, one can start from any initial condition and sample the binomial distributions in Eqs 26-28 together with the corresponding equations for the EZ nodes (overall 6 equations), over small time steps of size D ¼ � t=m t . Here, we used computer simulations to study the spread dynamics conditioned on seizure initiation at the EZ nodes at time t = 0. In addition, because of the separation of time-scales in the model, i.e.t r;i �t s;i , we considered different seizure spread events, separated by a long refractory period, as independent. Therefore, each realization starts with seizure initiation at the EZ nodes if EZ nodes are not inhibited (i.e. f(z 0 + E ez )>0), and the dynamics goes on until there is no seizure in the system.

Behavior of ψ − (apparent exponents)
For the behavior of ψ − (Eqs 29 and 31), the estimated exponentsb 0 andd 0 were negative. Negative values lead to a singularity in the function for ψ. Obviously, this type of singularity is not allowed. We present below a consistent formulation of the behavior of the order parameter spread size below the critical point w < w c where negative exponents were estimated. We follow the concept of apparent exponents as described in [32] in the context of estimating powerlaw probability density functions.
Letting x = w c − w, then we must have ψ(x = 0) = ψ c which in our model corresponds to a positive finite value. In addition, the response function w w ¼ @c @w ¼ À @c @x is expected to diverge at the critical point. Therefore, it follows that Based on the above conditions, consider next the following functional form for the order parameter: The function G must be in agreement with conditions 57 and 58. It is necessary that for finite positive A. In addition, it is necessary that y þb 0 > 0, which guarantees condition in Equation Eq 57. Using Eq 59, it is straightforward to show that the actual exponent governing the behavior of ψ near the critical point x ≳ 0 is b 0 ¼ y þb 0 .
The parameter x c in the above is a positive number that specifies the scale of x such that for x > x c , the function G approaches a constant value, thus leading to the observation of the apparent exponentb 0 . For small x < x c near x = 0, the function G scales as x θ , thus adding the value of θ to the apparent exponentb 0 .
From Eq 59 we can also calculate the response function whose behavior near the critical point x ≳ 0 is which diverges at the critical point if y þb < 1. Therefore, the intersection of conditions in Eqs 57 and 58 dictates 0 < y þb < 1: The above argumentation can be similarly applied to thed 0 exponent.